package com.vic.Report;

import java.util.Date;

public class CellData {
    /**
     * Initializes data. All input data must have the same length.
     */
    public CellData(Date[] rcv_time, double gen_vol[], double cell_vol[][], double gen_cur[]) {
        int len = rcv_time.length;
        if (gen_vol.length != len || gen_cur.length != len) {
            throw new IllegalArgumentException("All data must be of the same length");
        }
        for (int i = 0; i < cell_vol.length; ++i) {
            if (cell_vol[i].length != len) {
                throw new IllegalArgumentException("All data must be of the same length");
            }
        }
        this.rcv_time = rcv_time;
        this.gen_vol = gen_vol;
        this.cell_vol = cell_vol;
        this.gen_cur = gen_cur;
    }

    private final Date rcv_time[];
    private final double gen_vol[];
    private final double cell_vol[][];
    private final double gen_cur[];

    private static double sqr(double x) {
        return x * x;
    }

    private static double sigmoid(double x, double mu, double sigma) {
        double norm = (x - mu) / sigma;
        return 1. / (1 + Math.exp(-norm));
    }

    /**
     * Estimate the start and end point of discharge.
     * NOTE: the algorithm is not perfect, and will fail to extract the discharge interval on unusual data inputs.
     * The returned start point denotes the last sample before discharge begins,
     * and the end point denotes the last point before discharge ends.
     *
     * @return two-element array containing the start index and end index of discharge
     */
    public int[] estimateDischargeStartAndEnd() {
        double y[] = gen_vol;
        double i[] = gen_cur;
        if (y.length != i.length) {
            throw new IllegalArgumentException("Data must have the same length");
        }
        if (y.length <= 2) {
            int result[] = new int[2];
            result[0] = 0;
            result[1] = Math.max(0, y.length - 1);
            return result;
        }
        int len = y.length;
        double disc_cur_std = 0, delta_cur_std = 0, vol_mean = 0, vol_std = 0, delta_vol_std = 0, max_vol = 0;
        for (int k = 0; k < len; ++k) {
            if (i[k] < 0) {
                disc_cur_std += sqr(i[k]);
            }
            if (k > 0) {
                delta_cur_std += sqr(i[k] - i[k - 1]);
                delta_vol_std += sqr(y[k] - y[k - 1]);
                max_vol = Math.max(max_vol, y[k]);
            } else {
                max_vol = y[k];
            }
            vol_mean += y[k];
            vol_std += sqr(y[k]);
        }
        disc_cur_std = Math.sqrt(disc_cur_std / len);
        delta_cur_std = Math.sqrt(delta_cur_std / (len - 1));
        vol_mean /= len;
        vol_std = Math.sqrt(vol_std / len - sqr(vol_mean));
        delta_vol_std = Math.sqrt(delta_vol_std / (len - 1));
        double start_likelihood[] = new double[len];
        double end_likelihood[] = new double[len];
        for (int k = 0; k < len; ++k) {
            start_likelihood[k] = sigmoid(Math.min(0, i[k]), 0, disc_cur_std)                                           // i_k is positive
                * sigmoid(k > 0 ? Math.min(0, i[k] - i[k - 1]) : 0, 0, delta_cur_std * 2)                               // i_k - i_{k-1} is non-negative
                * sigmoid(k > 0 ? Math.min(0, y[k] - y[k - 1]) : 0, 0, delta_vol_std * 2)                               // y_k - y_{k-1} is non-negative
                * sigmoid(k < len - 1 ? i[k] - i[k + 1] : 0, 0, (delta_cur_std + disc_cur_std) * .5)                    // i_{k+1} - i_k is negative
                * sigmoid(k < len - 1 ? y[k] - y[k + 1] : 0, 0, (delta_vol_std + vol_std) * .5)                         // y_{k+1} - y_k is negative
                * sigmoid(y[k], vol_mean, vol_std)                                                                      // y_k is high
                * (sigmoid(k < len - 1 ? -Math.min(0, i[k + 1]) : -Math.min(0, i[k]), 0, disc_cur_std) - 0.5);          // i_{k+1} is negative
            end_likelihood[k] = (sigmoid(-Math.min(0, i[k]), 0, disc_cur_std) - 0.5)                                    // i_k is negative
                * sigmoid(k > 0 ? Math.min(0, i[k - 1] - i[k]) : 0, 0, delta_cur_std * 2)                               // i_k - i_{k-1} is non-positive
                * sigmoid(k > 0 ? Math.min(0, y[k - 1] - y[k]) : 0, 0, delta_vol_std * 2)                               // y_k - y_{k-1} is non-positive
                * sigmoid(k < len - 1 ? i[k + 1] - i[k] : -Math.min(0, i[k]), 0, (delta_cur_std + disc_cur_std) * .5)   // i_{k+1} - i_k is positive
                * sigmoid(k < len - 1 ? y[k + 1] - y[k] : max_vol - y[k], 0, (delta_vol_std + vol_std) * .5)            // y_{k+1} - y_k is positive
                * (1 - sigmoid(y[k], vol_mean, vol_std))                                                                // y_k is low
                * sigmoid(k < len - 1 ? Math.min(0, i[k + 1]) : 0, 0, disc_cur_std);                                    // i_{k+1} is positive
        }
        int max_end_likelihood_index = len - 1;
        double max_end_likelihood = end_likelihood[max_end_likelihood_index];
        int max_likelihood_start = len - 2;
        int max_likelihood_end = len - 1;
        double max_interval_likelihood = start_likelihood[max_likelihood_start] * end_likelihood[max_likelihood_end];
        for (int k = len - 3; k >= 0; --k) {
            if (end_likelihood[k + 1] > max_end_likelihood) {
                max_end_likelihood = end_likelihood[k + 1];
                max_end_likelihood_index = k + 1;
            }
            if (start_likelihood[k] * max_end_likelihood > max_interval_likelihood) {
                max_interval_likelihood = start_likelihood[k] * max_end_likelihood;
                max_likelihood_start = k;
                max_likelihood_end = max_end_likelihood_index;
            }
        }
        int result[] = new int[2];
        result[0] = max_likelihood_start;
        result[1] = max_likelihood_end;
        return result;
    }

    private static double CoulombCounting(Date t[], double i[], int start, int end) {
        double ampere_millisec = 0;
        for (int k = start; k < end; ++k) {
            ampere_millisec += -i[k] * (t[k + 1].getTime() - t[k].getTime());
        }
        return ampere_millisec / 3600000d;
    }

    /**
     * Estimate pack capacity and individual cell capacity using the given start and end point of discharge.
     * NOTE: for cells that did not discharge below the cutoff voltage, the capacity is estimated assuming a
     * third-order polynomial relationship between the final voltage and the cell capacity. This model is not
     * verified, so the precision of such estimation should not be relied upon.
     *
     * @param disc_start index of the last sample before discharge starts
     * @param disc_end index of the last sample before discharge ends
     *
     * @return a double array, containing the pack capacity as its first value, and the individual cell capacity
     * as the following values.
     */
    public double[] estimateCellCapacity(int disc_start, int disc_end) {
        int cell_cnt = get_cell_cnt();
        double capacity[] = new double[cell_cnt + 1];
        Date t[] = rcv_time;
        double i[] = gen_cur;
        int duration_start = Math.min(disc_start + 1, disc_end - 1);
        capacity[0] = CoulombCounting(t, i, duration_start, disc_end);
        double cutoff_voltage = 0, cutoff_volt_std = 0;
        for (int j = 0; j < cell_cnt; ++j) {
            double vol = get_cell_vol(j + 1)[disc_end];
            cutoff_voltage += vol;
            cutoff_volt_std += sqr(vol);
        }
        cutoff_voltage /= cell_cnt;
        cutoff_volt_std = Math.sqrt(cutoff_volt_std / cell_cnt - sqr(cutoff_voltage));
        int filtered_cnt = 0;
        double filtered_cutoff = 0;
        for (int j = 0; j < cell_cnt; ++j) {
            double vol = get_cell_vol(j + 1)[disc_end];
            if (vol >= cutoff_voltage - cutoff_volt_std * 3 && vol <= cutoff_voltage + cutoff_volt_std * 3) {
                filtered_cutoff += vol;
                ++filtered_cnt;
            }
        }
        if (filtered_cnt > 0) {
            filtered_cutoff /= filtered_cnt;
        } else {
            filtered_cutoff = cutoff_voltage;
        }
        double slant = 0.25;
        int known_cnt = 0;
        for (int j = 0; j < cell_cnt; ++j) {
            double yj[] = get_cell_vol(j + 1);
            if (yj[disc_end] <= filtered_cutoff) {
                int cutoff_ind = duration_start;
                for (int k = disc_end; k >= duration_start; --k) {
                    if (yj[k] >= filtered_cutoff * ((double) (k - duration_start) / (disc_end - duration_start) * slant - slant + 1)) {
                        cutoff_ind = k;
                        break;
                    }
                }
                capacity[j + 1] = CoulombCounting(t, i, duration_start, cutoff_ind);
                ++known_cnt;
            } else {
                capacity[j + 1] = Double.NaN;
            }
        }
        double x[];
        if (known_cnt < 2) {
            x = null;
        } else {
            // Assuming C_k = a*v_k^3 + b*v_k + c, identify the model using known final voltage and capacity values,
            // and then estimate the unknown capacity with the identified model
            double s1 = 0, s2 = 0, s3 = 0, s4 = 0, s6 = 0, sc0 = 0, sc1 = 0, sc3 = 0, vm1 = 0;
            for (int j = 0; j < cell_cnt; ++j) {
                double v1 = get_cell_vol(j + 1)[disc_end];
                double cap = capacity[j + 1];
                if (!Double.isNaN(cap)) {
                    double v2 = v1 * v1;
                    double v3 = v2 * v1;
                    s1 += v1;
                    s2 += v2;
                    s3 += v3;
                    s4 += v2 * v2;
                    s6 += v3 * v3;
                    sc0 += cap;
                    sc1 += v1 * cap;
                    sc3 += v3 * cap;
                }
                vm1 = Math.max(vm1, v1);
            }
            double H[], c[], A[], b[];
            if (known_cnt > 2) {
                H = new double[9];
                H[0] = s6;
                H[1] = H[3] = s4;
                H[2] = H[6] = s3;
                H[4] = s2;
                H[5] = H[7] = s1;
                H[8] = known_cnt;
                c = new double[3];
                c[0] = sc3;
                c[1] = sc1;
                c[2] = sc0;
                A = new double[15];
                double co1 = filtered_cutoff, co3 = co1 * co1 * co1;
                double vm3 = vm1 * vm1 * vm1;
                A[0] = A[6] = A[12] = A[13] = 1;
                A[1] = A[2] = A[5] = A[7] = A[10] = A[11] = 0;
                A[3] = co3;
                A[4] = -vm3;
                A[8] = co1;
                A[9] = -vm1;
                A[14] = -1;
                b = new double[5];
                b[0] = b[1] = b[2] = 0;
                b[3] = capacity[0];
                b[4] = capacity[0] * -2;
                x = new double[3];
                x[0] = 0;
                x[1] = 0;
                x[2] = capacity[0];
            } else {
                H = new double[4];
                H[0] = s2;
                H[1] = H[2] = s1;
                H[3] = known_cnt;
                c = new double[2];
                c[0] = sc1;
                c[1] = sc0;
                A = new double[8];
                double co1 = filtered_cutoff;
                A[0] = A[5] = A[6] = 1;
                A[1] = A[4] = 0;
                A[2] = co1;
                A[3] = -vm1;
                A[7] = -1;
                b = new double[4];
                b[0] = b[1] = 0;
                b[2] = capacity[0];
                b[3] = capacity[0] * -2;
                x = new double[2];
                x[0] = 0;
                x[1] = capacity[0];
            }
            ActiveSet solver = new ActiveSet();
            try {
                x = solver.solve(H, c, A, b, x);
            } catch (IllegalArgumentException ex) {
                x = null;
            }
        }
        if (x == null) {
            for (int j = 0; j < cell_cnt; ++j) {
                if (Double.isNaN(capacity[j + 1])) {
                    capacity[j + 1] = capacity[0];
                }
            }
        } else {
            for (int j = 0; j < cell_cnt; ++j) {
                if (Double.isNaN(capacity[j + 1])) {
                    double v1 = get_cell_vol(j + 1)[disc_end];
                    if (x.length == 3) {
                        capacity[j + 1] = (x[0] * v1 * v1 + x[1]) * v1 + x[2];
                    } else {
                        capacity[j + 1] = x[0] * v1 + x[1];
                    }
                }
            }
        }
        return capacity;
    }

    public Date[] get_rcv_time() {
        return rcv_time;
    }

    public double[] get_gen_vol() {
        return gen_vol;
    }

    public int get_cell_cnt() {
        return cell_vol.length;
    }

    public double[] get_cell_vol(int ind) {
        return cell_vol[ind - 1];
    }

    public double[] get_gen_cur() {
        return gen_cur;
    }
}
